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Many active fluid systems encountered in biology are set in total geometric confinement. Cyto- 
plasmic streaming in plant cells is a prominent and ubiquitous example, in which cargo-carrying 
molecular motors move along polymer filaments and generate coherent cell-scale flow. When fila- 
ments are not fixed to the cell periphery, a situation found both in vivo and in vitro, we observe 
that the basic dynamics of streaming are closely related to those of a dipolar microswimmer suspen- 
sion. This paradigm is used to demonstrate that confinement brings about a qualitative change in 
behavior; a linear stability analysis and numerical studies show that there is an activity threshold 
for spontaneous auto-circulation. Analysis of the long-time behavior reveals a phenomenon akin to 
defect separation in nematic liquid crystals, and a high-activity bifurcation to an oscillatory regime. 

PACS numbers: 87.16.Wd, 87.16.Ln, 47.63.-b, 47.54.-r 



Cytoplasmic streaming is the deliberate, driven mo- 
tion of the entire contents of large eukaryotic cells. It is 
effected by cargo-laden molecular motors walking along 
polymer filaments and entraining the surrounding fluid 
(Figure |TJl) ; the combined action of many of these mo- 
tors can generate flow speeds in excess of 100 /xm/s for 
certain freshwater algae. While inroads are being made 
into understanding its function [TJ [2] , surprisingly little 
is known about how it is initially established within cells. 

In a remarkable, yet apparently little-known investiga- 
tion into the development of streaming, Yotsuyanagi [3] 
in 1953 examined isolated droplets of cytoplasm forcibly 
extracted from algal cells. He observed a progression 
from isolated Brownian fluctuations to a coherent, global 
circulation of the entire droplet contents (Figure [TJd). 
However, we need not limit ourselves to ex vivo exper- 
iments: Kamiya [4 describes a similar blooming of ro- 
tational cyclosis in the development of Lilium pollen 
cells, and Jarosch [5] quantitatively analyzed the same 
disorder-to-order transition occurring within Allium cells 
over the course of a few hours. Based on these observa- 
tions, one is led to ask: is it possible that a simple self- 
organization process could lie at the heart of streaming? 

When the filaments are not locked in position, as is 
likely in Yotsuyanagi' s experiments, a cargo-carrying mo- 
tor walking on a free filament constitutes a force dipole. 
Therefore, these cytoplasmic dynamics belong to the bur- 
geoning field of active fluids. At its simplest, an active 
fluid is a continuum suspension of force dipoles interact- 
ing via short- and long-range forces, leading to a system 
like a liquid crystal, but with continuous injection of en- 
ergy at the microscale. Such systems can exhibit complex 
patterns and flows [6 , including asters and vortices [Ti- 
ll)] , laning [10] and density waves [9j [10] . While various 
complex short-range interactions can be included in these 
formulations, it is hydrodynamics that drives many of the 
pattern formation behaviors. 

Despite the ubiquity of relevant situations, of which 



streaming is a major example, the influence of total 
confinement is relatively little-studied. Recent work by 
Fiirthauer et al. [11 used spontaneous flow of generic 
active polar fluids to construct theoretically a 'Taylor- 
Couette motor' that can be used to do work against two 
cylindrical boundaries. However, this is topologically dis- 
tinct from a single confined chamber. Schaller et al. [12] 
underline the critical importance of long-range hydrody- 
namics in confined systems: swirling patterns were ob- 
served experimentally in a totally confined actin motility 
assay, but no such circulation was reproduced in sim- 
ple actin-substrate cellular automaton computer simula- 
tions. They concluded that confined flows are responsible 
for the formation and stability of the global circulation. 

Through theory and simulation, we demonstrate here 
that only two key ingredients are required to capture the 
spontaneous emergence of self-organized stable rotational 
flow in biological systems: activity and confinement. Our 
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FIG. 1. (color online). Cytoplasmic streaming in vivo and ex 
vivo, (a) A molecular motor attached to a vesicle (i) encoun- 
ters a filament, (ii) binds and walks along it, entraining fluid, 
before (iii) unbinding stochastically (b) A drop of cytoplasm 
extracted from a plant cell transitions from random Brownian 
fluctuations to ordered circulation [3]. 
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model assumes that short, rigid filaments are suspended 
in a Newtonian, zero Reynolds number fluid. The fila- 
ments are assumed to exert extensile, or 'pusher', dipo- 
lar forces on the fluid; this can be viewed as the effect 
of processive molecular motors landing randomly along 
a filament and walking toward one end, implying an av- 
erage motor location forward of the filament midpoint. 
Additionally, the suspension is taken to be dilute, so fil- 
aments interact via hydrodynamics only, and is confined 
within a no-slip sphere of diameter L. 

Working in d dimensions we generalize the standard 
kinetic approach to these systems [13] . The spatial and 
angular distribution function \P(x, p,£) of the filaments, 
where |p| = 1, satisfies a Smoluchowski equation 

^ = -V x .(x*)-V p .(p*) (1) 

where V x = <9/<9x and V p = (J— pp)-<9/<9p. The spatial 
and rotational fluxes are 

x = u + Vp - D (s) • V x log ^, 

p ee (J - pp) ( 7 E + W) p - D^Vp log *, 

where V is a self-advection speed, 7 G [—1,1] is a shape 
parameter (7 — >• 1 for a slender rod), is a spatial dif- 
fusion tensor and is a rotational diffusion constant. 
The fluid has velocity field u, rate-of-strain tensor E ee 
(Vu + Vu T )/2 and vorticity tensor W ee (Vu- Vu t )/2. 
The filament pusher stresslet of strength <j > generates 
a stress tensor X) ee — <j J p dp (pp — I/d)^ that drives 
fluid flow by the Stokes equation — /iV 2 u + VII = V • XI 
with viscosity \i and pressure II, subject to incompress- 
ibility V • u = 0. Confinement induces the no-slip 
boundary condition u = 0on|x|=L/2. 

While simulations of the full system are possible 
[T3HT5] . here we develop evolution equations for the pri- 
mary orientation moments [T6HT8] . Given the orienta- 
tional average ((/>) ee J p dp(/>\I/, define the concentration 
c ee (1), polar moment P ee (p) and nematic moment 
Q ee (pp — I/d). Equations of motion for these fields in 
terms of higher moments can then be derived by taking 
appropriate weighted integrals of Eq. [19] . 

We pare down complications by specializing to two di- 
mensions (d = 2), rodlike particles (7 = 1) and isotropic 
diffusion (D^ = D^I), and neglect self-advection (V = 
0). This last assumption decouples the c dynamics into 
pure advection-diffusion and eliminates all polar interac- 
tions, so we take a constant concentration c = cq and 
neglect P. However, the remaining Q dynamics still de- 
pends on the fourth moment contraction (pppp) : E, and 
a closure is needed. Typically this is done by taking the 
distribution ^ to be a functional purely of the first three 
moments, yielding a closure linear in Q [18 . In dense 
active systems this is permissible, owing to the presence 
of local interaction terms; here, however, it is the above 



fourth moment term which provides all stabilizing non- 
linearities, so greater care must be taken. Instead we 
adapt a closure of Hinch and Leal [20 to d = 2, yielding 

(pppp) :E^[4Q-E-Q + 2c(E Q + Q E) 

+ c 2 E -2IQ 2 : E] . 

After non-dimensionalizing by rescaling x — >• Lx, t — > 
(c L 2 //i)t, u (/i/coi)u, n -> (/i 2 /c L 2 )n, s 
(coL 2 1 /i 2 )5] and Q — >• coQ, the final model reads 

^ = d^V 2 Q - 4d^Q + \E - 2Q E Q (2) 

where D/Dt = d/dt + u • V, with non-dimensional 
diffusion constants d^ = (co//j)D^ and d^ = 
{cqL 2 I ' ii)D^ r \ This is subject to the Stokes equation 
— V 2 u + VII = — do V • Q and incompressibility V • u = 
with non-dimensional dipole stress ao = (c^L / ' ji) 2 a . The 
fluid boundary condition reads u = on |x| = 1/2. 
Among the variety of admissible boundary conditions on 
Q we focus here on the natural condition N • VQ = 0, 
where N is the boundary normal vector. Qualitatively 
similar results are found with fixed boundary-parallel or 
boundary-perpendicular conditions [T9] , 

The model ([2| has the structure of a Landau theory 
for the order parameter Q. As E is linear in the veloc- 
ity u, and u is (nonlocally) linear in <JoQ via the Stokes 
equation, the term (1/2) E oc 00 Q. It follows in the usual 
manner that there is an effective linear term in Q that will 
become positive for sufficiently large activity ctq relative 
to — 4d( r \ If this is sufficient to overcome the diffusive 
stabilization d^ then the amplitude of the ensuing insta- 
bility will be limited by the nonlinear term 2Q-E-Q oc Q 3 . 

We first seek a steady non-flowing axisymmetric state 
Q°. In polar coordinates a = (r,6) the tensor Laplacian 
of Q has primary components 

fv*n\ -rn - 1 9 ( d 9lA 4 ^ 
r or \ or J r z 

while the others follow from symmetry and the trace- 
lessness of Q. Eq. ^ therefore implies Q® r and Q® e 
each satisfy a (modified) Bessel equation in z = 2Ar, 
viz. z 2 d 2 Q° ra + zd z Q° ra - (z 2 + 4)Q° ra = 0, where 
A 2 ee dW/dW. Thus Q° ra ex J 2 (2Ar); since J 2 is mono- 
tonic, the boundary conditions imply Q° = everywhere. 

Now, perturb axisymmetrically: let Q = eR, e <C 1, 
and write u = evO, E = ee for the induced flow (which 
has no radial component by incompressibility). Seek an 
exponentially growing state such that dtR = sR. Then to 
O(e), the perturbation obeys sR = d^V 2 R-4d^R+±e. 
To determine e we employ the technique of Kruse et al. 
[7] and write the Stokes equation as V • (—11/ + 2e — 
ctqR) = V • X) tot = 0. The r-component determines II. 
The (9-component reads S r S^ t + (2/r)S^ t = 0, so for Yffi 
analytic at r = we find T i t r ° e t = 0, i.e. e r Q = (ao/2)R r Q. 
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FIG. 2. (color online). Numerical results beyond the sponta- 
neous circulation threshold. (a,b,c) Simulated schlieren tex- 
tures of nematic order director n (i.e. density plot of (n x n y ) 2 ). 
Lighter corresponds to diagonally-oriented filaments, darker 
to horizontal or vertical, (a) Steady circulation with a cen- 
tral spiral defect at low activity, (b) steady central defect 
separation into a pair of hyperbolic defects, (c) snapshot of 
oscillatory behavior with widely separated mobile defects, (d) 
Flow streamlines for low activity, showing circulation about 
the system center; darker streamlines indicate faster flow, (e) 
Enlargement of nematic director field structure in texture (b) , 
showing two hyperbolic defects, (f ) Flow streamlines for high- 
activity oscillation associated with texture (c) exhibiting off- 
center flow circulation. In all cases, = = 0.025. 

Finally, e rr — as there is no radial velocity component. 
The perturbation therefore satisfies 

d^CR rr = (4d (r) + s)R rr , (3) 
d^CR r e = (4d« + s - ^) Rre, (4) 

which are still of Bessel form. When s > —4d^ r \ Eq. ^ 
has a solution in terms of I2 , so boundary conditions im- 
ply R rr = 0. Now, let A = (4d^ + s - cr /4)/^ s ) and 
write Eq. Q as CR r e = \R r Q- For A > this again 
gives solutions in terms of I2 and so R r Q = 0. How- 
ever, for A < (i.e. <Jo sufficiently large) the solution is 
instead R r e oc J2W— Ar). Applying the boundary con- 



dition R' re (l/2) — yields the eigenvalue A = — 4?/q in 
terms of yo « 3.054, the first positive point satisfying 
J^il/o) = 0. This implies that the homogeneous disor- 
dered state is unstable to a spontaneously flowing mode 
when <j > a*, where (in physical units) 

a criterion we have verified numerically by full simula- 
tions of Eq. (2). To lend perspective, we consider typical 
values of the material properties. The stress amplitude 
can be expressed as a = f£, where / is the (typically pN) 
force exerted by motors and £ is the (typically jam) sepa- 
ration of the opposing forces of the stresslet. For micron- 
size rods we expect ~ 0.01 s _1 and ~ 10 -9 
cm 2 /s, so for system sizes L > 10 fim rotational diffusion 
dominates the parenthetical term in Eq. ([2]). Then for a 
fluid of the viscosity of water the instability will set in at 
concentrations greater than ~ 10 8 cm -3 , corresponding 
to a volume fraction well below 10 -3 . 

In numerical studies of the fully nonlinear dynamics 
we vary the dipolar activity ao while fixing the diffusion 
constants at d^ = d^ — 0.025, and use the eigende- 
composition Q = %n — //2), where the order param- 
eter S and (headless) director n are the degree of local 
alignment and the average alignment direction, respec- 
tively. For sufficiently weak activity above <r*, a stable 
steady state emerges of circulation about the system cen- 
ter (Figure |2|i&d). The spiral pattern of the nematic 
director field is reminiscent of the predictions of Kruse 
et al. [7 for unconfined active systems [see also [HI [21] . 
As ao is increased, stronger contributions emerge from 
higher radial modes in the spectrum of the order param- 
eter. Indeed, expanding the order parameter S(r, 9) in a 
Fourier series as S(r, 6) = S' n (r)e m6> , the axisymmet- 
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FIG. 3. (color online). Details of the bifurcation to circu- 
lation, (a) Numerically-evaluated steady state amplitudes 
I Sq 171 ^ I of Bessel series expansion for axisymmetric part of or- 
der parameter S at varying activity ao, with d^ — d^ = 
0.025. (b-d) Profiles of S at indicated points (b-d) in (a). 
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FIG. 4. (color online). Secondary bifurcation to oscilla- 
tory dynamics, (a) Amplitude of oscillation of the veloc- 
ity autocorrelation function H(r) as a function of <to, with 
c/( r ) = = 0.025, showing a bifurcation from steady defect 
separation to oscillatory behavior at a critical value of cr . (b) 
S(r) for do = 32 showing periodic oscillatory behavior, (c) 
Position of flow circulation center over time for ao = 32.75 
during one oscillation period. 

ric stability analysis suggests an n = mode expansion of 
the form S (r) Em=o S^ m) J 2 (2y m r) where J^Vm) 
and y m < ym+i- Mode amplitudes then read 

^o m) = T7 - d0 drrJ 2 (2y m r)S(r,0) 

JVm Jo JO 

with normalization M m = (7r/4)(l-4/^)[J 2 (?/ m )] 2 . Fig- 
ure [3] shows the steady-state values of the first three mode 
amplitudes \Sq°^ |, ISq 1 ^ |, | ^q 2 ^ | as functions of cr . Observe 
the non-zero contributions from modes with m > in the 
steady state despite only the m = mode being initially 
excited when a /16 < yfd^ + 

At larger values of ctq, the steady state exhibits Re- 
ject separation: the central axisymmetric spiral defect in 
the nematic director field (with topological charge +1) 
splits into two closely spaced hyperbolic defects (each of 
charge +1/2). The system still possesses fluid circulation 
about the central axis, due to the symmetric positioning 
of the defects. Such a configuration is illustrated in Fig- 
ure [2}}&e. The emergence of the defect separation phe- 
nomenon is perhaps unsurprising if we make contact with 
classical liquid crystal theory; for approximately isolated 
defects, the free energy penalty per defect is proportional 
to the square of its topological charge [22] , rendering two 
+1/2 defects favorable over a single +1 spiral. Indeed, de 
las Heras et al. [23] recently investigated the equivalent 
confined setup for a microscopic two-dimensional liquid 
crystal and always encountered defect separation. 

As <Jo is increased beyond a new critical value, the sys- 
tem bifurcates into a regime of periodic oscillation, where 
the time symmetry has been broken and a steady state is 



now unstable. The +1/2 defect pair (Figure]^) execute 
periodic 'orbits' around each other, with the flow circu- 
lation center offset from the origin (Figure |2f ) and fol- 
lowing an approximately circular trajectory (Figure^). 
These states can be analyzed by examining the velocity 
autocorrelation function [12] 

= / (V(x,t)-V(x,t-T)) x \ 

" l j -\ <v(x,t).v(x,t)) x A 

where the temporal average is taken over late times when 
the oscillatory state is fully established. Extracting the 
amplitude A of oscillation of S (Figure [4]o) we numeri- 
cally determine a bifurcation diagram as a function of ctq 
as in Figure [4^i. There is a clear threshold for the onset 
of periodic oscillations. 

Motivated by principles of cytoplasmic streaming, we 
have constructed a clean, simple model for a dilute sus- 
pension of extensile force-generating filaments in total 
geometric confinement, and have demonstrated that the 
inclusion of elementary hydrodynamics is entirely suffi- 
cient to yield spontaneous self-organization behavior, in 
spite of the absence of more complex local interaction 
terms. In an experimental realization, the prediction of 
a critical activity for transition from quiescence to cir- 
culation can be tested by varying the chamber size or 
the motor activity, perhaps through temperature or ATP 
concentration. Modern realizations of the experiment of 
Yotsuyanagi will likely provide a wealth of information 
on this type of bifurcation. 
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